A 3D finite element beam model via COMSOL

In this example, we demonstrate non-intrusive computation of SSMs via a 3D finite element beam model from COMSOL.

Finding a 2D SSM for a 3D finite element beam
generate model with a given point load
Dynamical system setup
Linear Modal analysis and SSM setup

Finding a 2D SSM for a 3D finite element beam

This is an example of how to reconstruct a slow 2D SSM of a mechanical system using all phase space variables measurements. In this example, we consider the clamped-clamped 3D beam that is modelled in COMSOL Multiphysics.

clear all
close all

% First we need to establish the connection between COMSOL and MATLAB from
% COMSOL server. See details in: https://doc.comsol.com/5.6/doc/com.comsol.help.llmatlab/llmatlab_ug_start.5.04.html
run ../../../install.m

generate model with a given point load

f0 = 1e5;

[model, M, C, K, f_0, Null, Nullf, ud, outdof, u0, disp_idx] = build_model(f0);

Construction of M,C,K of the finite element model in COMSOL. We use the function provided in COMSOL to extract the internal stiffness force.

[fint,dfint] = get_fint(K,model,Null, Nullf, ud, u0, disp_idx);

n = length(M);
disp(['Number of degrees of freedom = ' num2str(n)])
disp(['Phase space dimensionality = ' num2str(2*n)])

Number of degrees of freedom = 7875
Phase space dimensionality = 15750

Dynamical system setup

We consider the forced system

$\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}+\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})=\epsilon\mathbf{f}^{ext}(\mathbf{\Omega}t),$

which can be written in the first-order form as

$\mathbf{B}\dot{\mathbf{z}}=\mathbf{Az}+\mathbf{F}(\mathbf{z})+\epsilon\mathbf{F}^{ext}(\mathbf{\phi}) \dot{\mathbf{\phi}}=\mathbf{\Omega}$

where

$\mathbf{z}=\left[\begin{array}{c}\mathbf{x}\\\dot{\mathbf{x}}\end{array}\right],\quad\mathbf{A}=\left[\begin{array}{cc}-\mathbf{K} & \mathbf{0}\\\mathbf{0} & \mathbf{M}\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}\mathbf{C} & \mathbf{M}\\\mathbf{M} & \mathbf{0}\end{array}\right],\quad\quad\mathbf{F}(\mathbf{z})=\left[\begin{array}{c}\mathbf{-\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})}\\\mathbf{0}\end{array}\right],\quad\mathbf{F}^{ext}(\mathbf{z},\mathbf{\phi})=\left[\begin{array}{c}\mathbf{f}^{ext}(\mathbf{\phi})\\\mathbf{0}\end{array}\right]$ .

DSorder = 2;
DS = DynamicalSystem(DSorder);
set(DS,'M',M,'C',C,'K',K);
set(DS.Options, 'Intrusion', 'none')

set(DS,'fnl_non',fint,'dfnl_non',dfint);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')
set(DS.Options,'outDOF',outdof);

We assume periodic forcing of the form

$\mathbf{f}^{ext}(\phi) = \mathbf{f}_0\cos(\phi)=\frac{\mathbf{f}_0}{2}e^{i\phi} + \frac{\mathbf{f}_0}{2}e^{-i\phi}$

Fourier coefficients of Forcing

epsilon = 0.01;
kappas = [-1; 1];
coeffs = [f_0 f_0]/2;
DS.add_forcing(coeffs, kappas,epsilon);

Linear Modal analysis and SSM setup

[V,D,W] = DS.linear_spectral_analysis();

Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 2.021039e-03
modal damping ratio for 2 mode is 1.970247e-03
modal damping ratio for 3 mode is 1.996889e-03
modal damping ratio for 4 mode is 3.018811e-03
modal damping ratio for 5 mode is 3.036779e-03
the left eigenvectors may be incorrect in case of asymmetry of matrices

 The first 10 nonzero eigenvalues are given as 
   1.0e+03 *

  -0.0007 + 0.3682i
  -0.0007 - 0.3682i
  -0.0020 + 1.0056i
  -0.0020 - 1.0056i
  -0.0021 + 1.0357i
  -0.0021 - 1.0357i
  -0.0058 + 1.9344i
  -0.0058 - 1.9344i
  -0.0059 + 1.9486i
  -0.0059 - 1.9486i

Choose Master subspace (perform resonance analysis)

S = SSM(DS);
set(S.Options, 'reltol', 1,'IRtol',0.02,'notation', 'multiindex','contribNonAuto',false,'COMPtype','second')
set(S.FRCOptions, 'nt', 2^7, 'nRho', 300, 'nPar', 150, 'nPsi', 200, 'rhoScale', 5 )
set(S.FRCOptions, 'method','continuation ep' )
set(S.FRCOptions, 'outdof',outdof)

masterModes = [1,2];
S.choose_E(masterModes);

(near) outer resonance detected for the following combination of master eigenvalues
     2     0
     3     0
     3     1
     4     1
     4     2
     5     2
     0     2
     0     3
     1     3
     1     4
     2     4
     2     5
     2     0
     3     0
     3     1
     4     1
     4     2
     5     2
     0     2
     0     3
     1     3
     1     4
     2     4
     2     5
     5     0
     6     0
     6     1
     0     5
     0     6
     1     6
     5     0
     6     0
     6     1
     0     5
     0     6
     1     6

These are in resonance with the follwing eigenvalues of the slave subspace
   1.0e+03 *

  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0058 + 1.9344i
  -0.0058 + 1.9344i
  -0.0058 + 1.9344i
  -0.0058 - 1.9344i
  -0.0058 - 1.9344i
  -0.0058 - 1.9344i
  -0.0059 + 1.9486i
  -0.0059 + 1.9486i
  -0.0059 + 1.9486i
  -0.0059 - 1.9486i
  -0.0059 - 1.9486i
  -0.0059 - 1.9486i

sigma_out = 7
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     4     3
     1     2
     2     3
     3     4

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+02 *

  -0.0074 + 3.6815i
  -0.0074 + 3.6815i
  -0.0074 + 3.6815i
  -0.0074 - 3.6815i
  -0.0074 - 3.6815i
  -0.0074 - 3.6815i

sigma_in = 7

choose frequency range around the first natural frequency

omega0 = imag(S.E.spectrum(1));
omegaRange = omega0*[0.95 1.05];

extract forced response curve

order = 3;
FRC = S.extract_FRC('freq',omegaRange,order);

*****************************************
Calculating FRC using SSM with master subspace: [1  2]
(near) outer resonance detected for the following combination of master eigenvalues
     2     0
     3     0
     3     1
     4     1
     4     2
     5     2
     0     2
     0     3
     1     3
     1     4
     2     4
     2     5
     2     0
     3     0
     3     1
     4     1
     4     2
     5     2
     0     2
     0     3
     1     3
     1     4
     2     4
     2     5
     5     0
     6     0
     6     1
     0     5
     0     6
     1     6
     5     0
     6     0
     6     1
     0     5
     0     6
     1     6

These are in resonance with the follwing eigenvalues of the slave subspace
   1.0e+03 *

  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 + 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0020 - 1.0056i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 + 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0021 - 1.0357i
  -0.0058 + 1.9344i
  -0.0058 + 1.9344i
  -0.0058 + 1.9344i
  -0.0058 - 1.9344i
  -0.0058 - 1.9344i
  -0.0058 - 1.9344i
  -0.0059 + 1.9486i
  -0.0059 + 1.9486i
  -0.0059 + 1.9486i
  -0.0059 - 1.9486i
  -0.0059 - 1.9486i
  -0.0059 - 1.9486i

sigma_out = 7
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     4     3
     1     2
     2     3
     3     4

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+02 *

  -0.0074 + 3.6815i
  -0.0074 + 3.6815i
  -0.0074 + 3.6815i
  -0.0074 - 3.6815i
  -0.0074 - 3.6815i
  -0.0074 - 3.6815i

sigma_in = 7
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:00:31
Estimated memory usage at order  2 = 2.21E+01 MB
Manifold computation time at order 3 = 00:01:37
Estimated memory usage at order  3 = 2.26E+01 MB

 Run='freqSubint1Order3.ep': Continue equilibria along primary branch.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          1.92e-01  5.21e+02    0.0    0.0    0.0
   1   1  1.00e+00  1.69e-01  1.38e-02  5.21e+02    0.0    0.0    0.0
   2   1  1.00e+00  1.76e-03  5.02e-06  5.21e+02    0.0    0.0    0.0
   3   1  1.00e+00  5.20e-07  3.11e-13  5.21e+02    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:00   5.2071e+02      1  EP      3.6815e+02   5.1670e-02   5.6919e+00   1.0000e-02
   10  00:00:00   5.1660e+02      2          3.6524e+02   2.1643e-02   6.0480e+00   1.0000e-02
   20  00:00:00   5.1160e+02      3          3.6171e+02   1.0582e-02   6.1691e+00   1.0000e-02
   30  00:00:00   5.0661e+02      4          3.5817e+02   6.8951e-03   6.2089e+00   1.0000e-02
   40  00:00:01   5.0161e+02      5          3.5463e+02   5.1038e-03   6.2282e+00   1.0000e-02
   50  00:00:01   4.9661e+02      6          3.5110e+02   4.0494e-03   6.2396e+00   1.0000e-02
   54  00:00:01   4.9469e+02      7  EP      3.4974e+02   3.7519e-03   6.2428e+00   1.0000e-02

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:01   5.2071e+02      8  EP      3.6815e+02   5.1670e-02   5.6919e+00   1.0000e-02
   10  00:00:01   5.2474e+02      9          3.7101e+02   8.6793e-02   5.0588e+00   1.0000e-02
   20  00:00:01   5.2574e+02     10          3.7173e+02   9.1861e-02   4.6271e+00   1.0000e-02
   22  00:00:01   5.2574e+02     11  SN      3.7173e+02   9.1644e-02   4.6026e+00   1.0000e-02
   22  00:00:01   5.2574e+02     12  FP      3.7173e+02   9.1644e-02   4.6026e+00   1.0000e-02
   30  00:00:01   5.2569e+02     13          3.7169e+02   8.9987e-02   4.4916e+00   1.0000e-02
   40  00:00:02   5.2435e+02     14          3.7075e+02   6.7502e-02   3.9592e+00   1.0000e-02
   48  00:00:02   5.2376e+02     15  SN      3.7034e+02   4.5721e-02   3.6571e+00   1.0000e-02
   48  00:00:02   5.2376e+02     16  FP      3.7034e+02   4.5721e-02   3.6571e+00   1.0000e-02
   50  00:00:02   5.2376e+02     17          3.7034e+02   4.4368e-02   3.6403e+00   1.0000e-02
   60  00:00:02   5.2385e+02     18          3.7040e+02   3.8542e-02   3.5699e+00   1.0000e-02
   70  00:00:02   5.2586e+02     19          3.7183e+02   1.9202e-02   3.3498e+00   1.0000e-02
   80  00:00:02   5.3086e+02     20          3.7536e+02   9.5891e-03   3.2450e+00   1.0000e-02
   90  00:00:02   5.3586e+02     21          3.7890e+02   6.4293e-03   3.2108e+00   1.0000e-02
  100  00:00:03   5.4086e+02     22  EP      3.8243e+02   4.8382e-03   3.1937e+00   1.0000e-02
Time spent on mapping the FRC to physical coordinates maptime 7.8981 seconds. 
Total time spent on FRC computation upto O(3) = 00:02:21

A 3D finite element beam model via COMSOL

Contents

Finding a 2D SSM for a 3D finite element beam

generate model with a given point load

Dynamical system setup

Linear Modal analysis and SSM setup